*******		REPLICATION FILES
*******		Climate Variability and Irregular Migration to the European Union 
*******		Global Environmental Change 
*******		Fabien Cottier and Idean Salehyan
*******		Replication: Sensitivity analysis (S13) spatial heteroscedasticity and autocorrelation consistent (sHAC) (Tables A.40-42, Figure A.29-30)
*******		This version: April 2021

* note: stata packages required for running models
* - coefplot
* - crossfold
* - plottig
* - ols_spatial_HAC



* IMPORTANT: this file REQUIRES the ols_spatial_HAC function (Hsiang 2010). The README file contains the link to the repository storing the Stata ado files on Solomon Hsiang research group website at Berkeley. THESE MUST BE INSTALLED PRIOR TO RUNNING THIS SCRIPT.


clear all

set more off
set scheme plottig
cd "path/to/replication/directory"

* load data
use "data_quarter.dta", clear 

* load stata risk scripts
qui do scripts/script_rr_spei_spatialHAC_012021

* duplicates data
duplicates list cowc year quarter
duplicates tag cowc year quarter, gen(_dupl)
duplicates drop cowc year quarter, force

* set time series
sort cowc year quarter
gen quarter_int = real(regexs(1)) if regexm(quarter,"([0-9]+)")
gen quarterS=(year-2005)*4+quarter_int
order cowc nationality year quarter quarter_int quarterS
tsset cowc quarterS
 
* gen additional variables
encode continent, gen(contN)
recode contN (5=.)


* generate lag quartely migration variables
by cowc: gen nmigrq_exbalk_1tl = nmigrq_exbalk[_n-1]
by cowc: gen nmigrq_exbalk_2tl = nmigrq_exbalk[_n-2]
by cowc: gen nmigrq_exbalk_3tl = nmigrq_exbalk[_n-3]
by cowc: gen nmigrq_exbalk_4tl = nmigrq_exbalk[_n-4]
by cowc: gen nmigrq_exbalk_ln_1tl = nmigrq_exbalk_ln[_n-1]
by cowc: gen nmigrq_exbalk_ln_2tl = nmigrq_exbalk_ln[_n-2]
by cowc: gen nmigrq_exbalk_ln_3tl = nmigrq_exbalk_ln[_n-3]
by cowc: gen nmigrq_exbalk_ln_4tl = nmigrq_exbalk_ln[_n-4]


* gen dummies variables for sample splitting: low agriculture / high agriculture
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* i.year if migr100_exbalk==1, ///
 fe vce(cluster cowc)

sum agriLab if e(sample) & year==2010, d
scalar agriMed=r(p50)
gen agri_H=0 if e(sample) & year==2010
recode agri_H (0=1) if agriLab-agriMed > 0 & year==2010
by cowc: replace agri_H = agri_H[21] if missing(agri_H)

* generate dummies variable for extreme weather (cut-off 10 percentile / 90 percentile)
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* wmean_speiy /// 
 i.year i.quarter_int if migr100_exbalk==1, ///
 fe vce(cluster cowc)

centile (wmean_speiy) if e(sample), centile (10 90)
gen spei_drought = 0 if wmean_speiy!=.
recode spei_drought (0=1) if wmean_speiy <= r(c_1) 
centile (wmean_speiy) if e(sample), centile (10 90)
gen spei_hrain = 0 if wmean_speiy!=.
recode spei_hrain (0=1) if wmean_speiy >=  r(c_2) 


* create dummies for countries, years and quarter + constant

* constant
gen const = 1

* dummies for fixed effects
tabulate cowc, gen(cowc_dum)
tabulate year, gen(year_dum)
tabulate quarter_int, gen(quarter_dum)

* drop specific baseline dummies
drop year_dum6 // year == 2009
drop quarter_dum1 // quarter == 1

* ols_spatial_HAC does not accept interacted variables
gen wmean_speiy_sq = wmean_speiy^2


*******		  Macros

global SPEI_Y0 wmean_speiy

global SPEI_SQ_Y0 wmean_speiy wmean_speiy_sq
 
global SPEI_Y0Y2 L(0 4 8).wmean_speiy

global SPEI_SQ_Y0Y2 L(0 4 8).wmean_speiy L(0 4 8).wmean_speiy_sq


*******		 Table A.40


* NULL Model // No SPEI term 
qui ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 const ///
  cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1, ///
  lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000)  ///
  dropvar
qui ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* const ///
  cowc_dum* year_dum* quarter_dum* if e(sample), ///
  lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
  dropvar
est sto tabA40_m0
* list countries included in sample 
unique cowc if e(sample)
tabulate nationality contN if e(sample)


* Model 1
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 const /// 
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA40_m1
* F test spei parameters & aic
test wmean_speiy


* Model 2
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_SQ_Y0 const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA40_m2
* F test spei parameters & aic
test wmean_speiy wmean_speiy_sq


* Model 3
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0Y2 const /// 
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA40_m3
* F test spei parameters & aic
test L0.wmean_speiy L4.wmean_speiy L8.wmean_speiy


* Model 4
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_SQ_Y0Y2 const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA40_m4
* F test spei parameters & aic
test L0.wmean_speiy L4.wmean_speiy L8.wmean_speiy ///
  L0.wmean_speiy_sq L4.wmean_speiy_sq L8.wmean_speiy_sq




*******		 Table A.41


* NULL Model // agrarian countries
qui ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 const ///
  cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==1, ///
  lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
  dropvar
qui ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* const /// 
  cowc_dum* year_dum* quarter_dum* if e(sample), ///
  lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
  dropvar
est sto tabA41_m0H
* list of countries
unique cowc if e(sample)
tabulate nationality contN if e(sample)

* NULL Model // non-agrarian countries
qui ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 const ///
  cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==0, ///
  lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
  dropvar
qui ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* const ///
  cowc_dum* year_dum* quarter_dum* if e(sample), ///
  lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
  dropvar
est sto tabA41_m0L
* list of countries
unique cowc if e(sample)
tabulate nationality contN if e(sample)



* Model 5 // agrarian countries
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==1, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA41_m5
* F test spei parameters & aic
test wmean_speiy

* Model 6 // non-agrarian countries
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0 const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==0, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA41_m6
* F test spei parameters & aic
test wmean_speiy


* Model 7 // agrarian countries
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_SQ_Y0 const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==1, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA41_m7
* F test spei parameters & aic
test wmean_speiy wmean_speiy_sq

* Model 8 // non-agrarian countries
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_SQ_Y0 const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==0, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA41_m8
* F test spei parameters & aic
test wmean_speiy wmean_speiy_sq


* Model 9 // agrarian countries
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0Y2 const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==1, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA41_m9
unique cowc if e(sample)
* F test spei parameters & aic
test L0.wmean_speiy L4.wmean_speiy L8.wmean_speiy

* Model 10 // non-agrarian countries
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_Y0Y2 const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==0, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA41_m10
unique cowc if e(sample)
* F test spei parameters & aic
test L0.wmean_speiy L4.wmean_speiy L8.wmean_speiy


* Model 11 // agrarian countriese
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_SQ_Y0Y2 const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==1, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA41_m11
* F test spei parameters & aic
test L0.wmean_speiy L4.wmean_speiy L8.wmean_speiy ///
  L0.wmean_speiy_sq L4.wmean_speiy_sq L8.wmean_speiy_sq

* Model 12 // non-agrarian countries
* low reliance on agriculture
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* $SPEI_SQ_Y0Y2 const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==0, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA41_m12
* F test spei parameters & aic
test L0.wmean_speiy L4.wmean_speiy L8.wmean_speiy ///
  L0.wmean_speiy_sq L4.wmean_speiy_sq L8.wmean_speiy_sq



*******		 Table A.42



* Model 13 // agrarian countries
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* spei_drought spei_hrain const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==1, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA42_m13
* F test spei parameters & aic
test spei_drought spei_hrain

* Model 14 // non-agrarian countries
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* spei_drought spei_hrain const  ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==0, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA42_m14
* F test spei parameters & aic
test spei_drought spei_hrain


* Model 15 // agrarian countries
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* L(0 4 8).spei_drought L(0 4 8).spei_hrain const ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==1, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA42_m15
* F test spei parameters & aic
test L0.spei_drought L4.spei_drought L8.spei_drought L0.spei_hrain L4.spei_hrain L8.spei_hrain

* Model 16 // non-agrarian countries
ols_spatial_HAC nmigrq_exbalk_ln nmigrq_exbalk_ln_* L(0 4 8).spei_drought L(0 4 8).spei_hrain const  ///
	cowc_dum* year_dum* quarter_dum* if migr100_exbalk==1 & agri_H==0, /// 
	lat(lat) lon(lon) timevar(quarterS) panelvar(cowc) lag(2) dist(1000) ///
	dropvar
est sto tabA42_m16
* F test spei parameters & aic
test L0.spei_drought L4.spei_drought L8.spei_drought L0.spei_hrain L4.spei_hrain L8.spei_hrain



*******		 Relative risk figures




* Figure A.29 (Model 2, Table A.40)

qui rr_est tabA40_m2 "quadratic"


* Figure A.30 (Model 7 & 8, Table A.41)

qui rr_est_Sp tabA41_m7 tabA41_m8 "quadratic"

